Is canopy treatment or canopy cover a stronger predictor of the environmental conditions below the shrub canopy?
Notes: * response variable: The variable you want to test. (e.g.,. PAR, TDR, nncover) * explanatory variables: The variables that can explain the outcome. (e.g., treatment, mafacover, site) * data frame: Include all the data you’ll be working with. (e.g., combine par and mafa_cover_all) * formula: a formula is a two-sided linear formula object describing both the fixed-effects and random effects part of the model, with the response on the left of a ~ operator and the terms, separated by + operators, on the right. * random effects: Random-effects terms are distinguished by vertical bars (|) separating expressions for design matrices from grouping factors. Two vertical bars (||) can be used to specify multiple uncorrelated random effects for the same grouping variable. + Because of the way it is implemented, the ||-syntax works only for design matrices containing numeric (continuous) predictors; to fit models with independent categorical effects, see dummy or the lmer_alt function from the afex package.)
lme.parmod <- lmer(normalized_par ~ treatment + mafacover + site + (1|date) +(1|plot_number),
data = parmod)
MOD <- lme.parmod
# Calculate model's residual and fitted values
F1 <- fitted(MOD)
E1 <- residuals(MOD, "pearson")
# See residuals as a histogram, look for normal distribution
hist(E1)
# Graph residuals against a normal distribution, look for all points to be on 1:1 line
#need to call next two lines together
qqnorm(E1, xlab = "Theoretical Quantiles", ylab = "Pearson Residuals")
qqline(E1)
# Graph residuals vs fitted values to inspect homogeneity of variance, look for a random scattering of plots
# run together 2 lines
plot(x=F1, y=E1, xlab="Fitted Model Values", ylab="Pearson Residuals", main = "Variance")
abline(h = 0) # looking for random scatter
# Check for crazy outliers, look for lines MUCH taller than the rest of data
plot(cooks.distance(MOD))
# Plot a box plot for each factor, look for equal variance amount levels along the zero line
# Box plots need to be assessed for every fixed effect (categorical factor) for your model.
# Random effects need the qqline just like normality.
# So, only site and treatment get box plots.
# MAFA cover is a covariate (continuous variable).
boxplot(E1 ~ treatment, data = parmod)
abline(h = 0)
boxplot(E1 ~ site, data = parmod)
abline(h = 0)
# Plot each random effect in the model - things in the (1|_), look for dots along 1:1 line
qqnorm(ranef(MOD)$date[,1], main = "Plot Residuals")
qqline(ranef(MOD)$date[,1]) # unique dates
qqnorm(ranef(MOD)$plot_number[,1], main = "Plot Residuals")
qqline(ranef(MOD)$plot_number[,1]) # unique plots
3. statistics
# print
print(lme.parmod)
## Linear mixed model fit by REML ['lmerModLmerTest']
## Formula: normalized_par ~ treatment + mafacover + site + (1 | date) +
## (1 | plot_number)
## Data: parmod
## REML criterion at convergence: -247.8414
## Random effects:
## Groups Name Std.Dev.
## plot_number (Intercept) 0.09588
## date (Intercept) 0.08129
## Residual 0.13089
## Number of obs: 288, groups: plot_number, 36; date, 8
## Fixed Effects:
## (Intercept) treatmenthalf treatmentno mafacover site
## 0.413527 -0.004588 -0.011217 -0.003168 0.035283
summary(lme.parmod) # p-values
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: normalized_par ~ treatment + mafacover + site + (1 | date) +
## (1 | plot_number)
## Data: parmod
##
## REML criterion at convergence: -247.8
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.9319 -0.7023 -0.0626 0.6010 3.1787
##
## Random effects:
## Groups Name Variance Std.Dev.
## plot_number (Intercept) 0.009192 0.09588
## date (Intercept) 0.006608 0.08129
## Residual 0.017131 0.13089
## Number of obs: 288, groups: plot_number, 36; date, 8
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 4.135e-01 8.065e-02 4.954e+01 5.128 4.87e-06 ***
## treatmenthalf -4.588e-03 4.610e-02 3.458e+01 -0.100 0.92129
## treatmentno -1.122e-02 4.892e-02 3.856e+01 -0.229 0.81987
## mafacover -3.168e-03 9.918e-04 1.116e+02 -3.194 0.00183 **
## site 3.528e-02 4.077e-02 3.995e+01 0.866 0.39193
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) trtmnth trtmntn mafcvr
## treatmnthlf -0.074
## treatmentno 0.009 0.572
## mafacover -0.541 -0.333 -0.459
## site -0.841 -0.164 -0.226 0.492
anova(lme.parmod) # fixed level effects, pvalues, etc.
## Type III Analysis of Variance Table with Satterthwaite's method
## Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## treatment 0.000926 0.000463 2 34.385 0.0270 0.973361
## mafacover 0.174745 0.174745 1 111.567 10.2005 0.001825 **
## site 0.012833 0.012833 1 39.948 0.7491 0.391927
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# scatterplot: par x mafacover
ggplot(data = parmod, aes(x = mafacover, y = normalized_par)) +
geom_point() +
geom_smooth(method="lm")+
theme_bw()
## `geom_smooth()` using formula = 'y ~ x'
# scatterplot: par x site
ggplot(data = parmod, aes(x = site, y = normalized_par)) +
geom_boxplot(aes(fill = site)) +
theme_bw()
## Warning: Continuous x aesthetic
## ℹ did you forget `aes(group = ...)`?
## Warning: The following aesthetics were dropped during statistical transformation: fill
## ℹ This can happen when ggplot fails to infer the correct grouping structure in
## the data.
## ℹ Did you forget to specify a `group` aesthetic or to convert a numerical
## variable into a factor?
# scatterplot: mafacover x par x treatment
ggplot(data = parmod, aes(x = mafacover, y = normalized_par)) +
geom_point(aes(color = treatment)) +
theme_bw()
# boxplot: mafacover x site
ggplot(data = parmod, aes(x = site, y = normalized_par)) +
geom_boxplot(aes(fill = site)) +
theme_bw() +
theme(legend.position = "none")
## Warning: Continuous x aesthetic
## ℹ did you forget `aes(group = ...)`?
## The following aesthetics were dropped during statistical transformation: fill
## ℹ This can happen when ggplot fails to infer the correct grouping structure in
## the data.
## ℹ Did you forget to specify a `group` aesthetic or to convert a numerical
## variable into a factor?
parmod_treatment <- aov(normalized_par ~ treatment, data = parmod)
TukeyHSD(parmod_treatment)
## Tukey multiple comparisons of means
## 95% family-wise confidence level
##
## Fit: aov(formula = normalized_par ~ treatment, data = parmod)
##
## $treatment
## diff lwr upr p adj
## half-full -0.05365265 -0.11695323 0.009647926 0.1147974
## no-full -0.08296590 -0.14626648 -0.019665321 0.0062500
## no-half -0.02931325 -0.09261383 0.033987332 0.5204084
plot(TukeyHSD(parmod_treatment))
lme.tdrmod <- lmer(soil_moisture ~ treatment + mafacover + site + (1|date) +(1|plot_number),
data = tdrmod)
MOD <- lme.tdrmod
# Calculate model's residual and fitted values
F1 <- fitted(MOD)
E1 <- residuals(MOD, "pearson")
# See residuals as a histogram, look for normal distribution
hist(E1)
# Graph residuals against a normal distribution, look for all points to be on 1:1 line
#need to call next two lines together
qqnorm(E1, xlab = "Theoretical Quantiles", ylab = "Pearson Residuals")
qqline(E1)
# Graph residuals vs fitted values to inspect homogeneity of variance, look for a random scattering of plots
# run together 2 lines
plot(x=F1, y=E1, xlab="Fitted Model Values", ylab="Pearson Residuals", main = "Variance")
abline(h = 0) # looking for random scatter
# Check for crazy outliers, look for lines MUCH taller than the rest of data
plot(cooks.distance(MOD))
# Plot a box plot for each factor, look for equal variance amount levels along the zero line
# Box plots need to be assessed for every fixed effect (categorical factor) for your model.
# Random effects need the qqline just like normality.
# So, only site and treatment get box plots.
# MAFA cover is a covariate (continuous variable).
boxplot(E1 ~ treatment, data = tdrmod)
abline(h = 0)
boxplot(E1 ~ site, data = tdrmod)
abline(h = 0)
# Plot each random effect in the model, look for dots along 1:1 line
qqnorm(ranef(MOD)$date[,1], main = "Plot Residuals")
qqline(ranef(MOD)$date[,1])
qqnorm(ranef(MOD)$plot_number[,1], main = "Plot Residuals")
qqline(ranef(MOD)$plot_number[,1])
3. statistics
# print
print(lme.tdrmod)
## Linear mixed model fit by REML ['lmerModLmerTest']
## Formula: soil_moisture ~ treatment + mafacover + site + (1 | date) + (1 |
## plot_number)
## Data: tdrmod
## REML criterion at convergence: 522.5547
## Random effects:
## Groups Name Std.Dev.
## plot_number (Intercept) 1.372
## date (Intercept) 7.207
## Residual 2.249
## Number of obs: 108, groups: plot_number, 36; date, 4
## Fixed Effects:
## (Intercept) treatmenthalf treatmentnone mafacover site
## 21.81315 -0.86322 -3.15994 0.01808 -4.06263
summary(lme.tdrmod)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: soil_moisture ~ treatment + mafacover + site + (1 | date) + (1 |
## plot_number)
## Data: tdrmod
##
## REML criterion at convergence: 522.6
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.07912 -0.53157 0.00198 0.51062 2.23084
##
## Random effects:
## Groups Name Variance Std.Dev.
## plot_number (Intercept) 1.884 1.372
## date (Intercept) 51.934 7.207
## Residual 5.060 2.249
## Number of obs: 108, groups: plot_number, 36; date, 4
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 21.81315 3.90031 4.07414 5.593 0.00475 **
## treatmenthalf -0.86322 0.82965 35.78406 -1.040 0.30511
## treatmentnone -3.15994 0.92207 41.74189 -3.427 0.00138 **
## mafacover 0.01808 0.02545 79.56419 0.710 0.47961
## site -4.06263 0.83253 55.27825 -4.880 9.4e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) trtmnth trtmntn mafcvr
## treatmnthlf -0.006
## treatmentnn 0.045 0.591
## mafacover -0.233 -0.368 -0.548
## site -0.352 -0.199 -0.297 0.542
anova(lme.tdrmod) # A significant p-value indicates that
## Type III Analysis of Variance Table with Satterthwaite's method
## Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## treatment 66.942 33.471 2 36.100 6.6149 0.003568 **
## mafacover 2.553 2.553 1 79.564 0.5045 0.479615
## site 120.493 120.493 1 55.278 23.8131 9.403e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# one or more significant differences exist between group means.
tdrmod_site <- aov(soil_moisture ~ site*treatment, data = tdrmod)
anova(tdrmod_site)
## Analysis of Variance Table
##
## Response: soil_moisture
## Df Sum Sq Mean Sq F value Pr(>F)
## site 1 344.6 344.57 7.7568 0.00638 **
## treatment 2 154.9 77.45 1.7434 0.18009
## site:treatment 2 4.6 2.29 0.0516 0.94971
## Residuals 102 4531.0 44.42
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
TukeyHSD(tdrmod_site)
## Warning in replications(paste("~", xx), data = mf): non-factors ignored: site
## Warning in replications(paste("~", xx), data = mf): non-factors ignored: site,
## treatment
## Warning in TukeyHSD.aov(tdrmod_site): 'which' specified some non-factors which
## will be dropped
## Tukey multiple comparisons of means
## 95% family-wise confidence level
##
## Fit: aov(formula = soil_moisture ~ site * treatment, data = tdrmod)
##
## $treatment
## diff lwr upr p adj
## half-full -0.6462963 -4.382655 3.090062 0.9110153
## none-full -2.8011574 -6.537516 0.935201 0.1803794
## none-half -2.1548611 -5.891219 1.581497 0.3595155
plot(TukeyHSD(tdrmod_site))
## Warning in replications(paste("~", xx), data = mf): non-factors ignored: site
## Warning in replications(paste("~", xx), data = mf): non-factors ignored: site,
## treatment
## Warning in TukeyHSD.aov(tdrmod_site): 'which' specified some non-factors which
## will be dropped
4. figures
# boxplot: TDR by treatment by site
ggplot(data = tdrmod, aes(x = treatment, y = soil_moisture)) +
geom_boxplot(aes(fill = site)) +
theme_bw()
## Warning: The following aesthetics were dropped during statistical transformation: fill
## ℹ This can happen when ggplot fails to infer the correct grouping structure in
## the data.
## ℹ Did you forget to specify a `group` aesthetic or to convert a numerical
## variable into a factor?
# boxplot: TDR by site by treatment
ggplot(data = tdrmod, aes(x = site, y = soil_moisture)) +
geom_boxplot(aes(fill = treatment)) +
scale_fill_manual(values = c('#F6EECF', '#B09175', '#291611')) +
labs(title = "Soil Moisture by Site and Treatment", x = "Site", y = "PAR") +
theme_bw()
# scatterplot: tdr x mafacover
ggplot(data = tdrmod, aes(x = mafacover, y = soil_moisture)) +
geom_point() +
geom_smooth(method="lm")+
theme_bw()
## `geom_smooth()` using formula = 'y ~ x'
#need to convert df to a tibble
#as_tibble(nnmod)
#print(nn)
#names(nnmod)
#print(nnmod)
#str(nnmod)
lme.nnmod <- lmer(nncover ~ treatment + mafacover + site + (1|date) + (1|plot_number),
data = nnmod)
MOD <- lme.nnmod
# Calculate model's residual and fitted values
F1 <- fitted(MOD)
E1 <- residuals(MOD, "pearson")
# See residuals as a histogram, look for normal distribution
hist(E1)
# Graph residuals against a normal distribution, look for all points to be on 1:1 line
#need to call next two lines together
qqnorm(E1, xlab = "Theoretical Quantiles", ylab = "Pearson Residuals")
qqline(E1)
# Graph residuals vs fitted values to inspect homogeneity of variance, look for a random scattering of plots
# run together 2 lines
plot(x=F1, y=E1, xlab="Fitted Model Values", ylab="Pearson Residuals", main = "Variance")
abline(h = 0) # looking for random scatter
# Check for crazy outliers, look for lines MUCH taller than the rest of data
plot(cooks.distance(MOD))
# Plot a box plot for each factor, look for equal variance amount levels along the zero line
# Box plots need to be assessed for every fixed effect (categorical factor) for your model.
# Random effects need the qqline just like normality.
# So, only site and treatment get box plots.
# MAFA cover is a covariate (continuous variable).
boxplot(E1 ~ treatment, data = nnmod)
abline(h = 0)
boxplot(E1 ~ site, data = nnmod)
abline(h = 0)
# Plot each random effect in the model, look for dots along 1:1 line
qqnorm(ranef(MOD)$date[,1], main = "Plot Residuals")
qqline(ranef(MOD)$date[,1])
qqnorm(ranef(MOD)$plot_number[,1], main = "Plot Residuals")
qqline(ranef(MOD)$plot_number[,1])
3. statistics
# boxplot: nncover x treatment x site
ggplot(data = nnmod, aes(x = treatment, y = nncover)) +
geom_boxplot(aes(fill = site)) +
geom_smooth(method="lm") +
labs(title = "Non-native cover by Treatment and Site", x = "Treatment", y = "Non-native cover (%)") +
scale_fill_manual(values = c( "#DCC27A", "#F19B34")) +
# scale_fill_manual(values = c('#D3E3CA', '#92A587')) +
theme_bw()
## Warning: The following aesthetics were dropped during statistical transformation: fill
## ℹ This can happen when ggplot fails to infer the correct grouping structure in
## the data.
## ℹ Did you forget to specify a `group` aesthetic or to convert a numerical
## variable into a factor?
## `geom_smooth()` using formula = 'y ~ x'
# boxplot: mafacover x site x treatment
ggplot(data = nnmod, aes(x = site, y = nncover)) +
geom_boxplot(aes(fill = treatment)) +
# scale_fill_manual(values = c('#F19B34', '#9F7E75', '#92A587')) +
scale_fill_manual(values = c('#F6EECF', '#B09175', '#291611')) +
labs(title = "Non-native cover by Site and Treatment", x = "Site", y = "Non-native cover (%)") +
geom_smooth(method="lm") +
theme_bw()
## `geom_smooth()` using formula = 'y ~ x'
# scatterplot: nncover x mafacover
ggplot(data = nnmod, aes(x = mafacover, y = nncover)) +
geom_point() +
geom_smooth(method="lm")+
theme_bw()
## `geom_smooth()` using formula = 'y ~ x'
#Question 2: How does seedling survival change by PAR, TDR, and mafacover? Is canopy treatment or canopy cover a stronger predictor of the environmental conditions below the shrub canopy?
df exploration * names(df) will give you the column names in the df
statistics * lmer: Fit Linear Mixed-Effects Models + Fit a linear mixed-effects model (LMM) to data, via REML or maximum likelihood.
end of code